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Abstract 

A mixed classical-quantum density functional theory is used to calculate pair 
correlations and the free energy of a spin-polarized Hydrogen plasma. A transition 
to an atomic insulator phase is estimated to occur around r s = 2.5 at T = lO 4 -^, 
and a pressure P « 0.5Mbar. Spin polarization is imposed to prevent the formation 
of Hi molecules. 



2 



Although Hydrogen is generally considered to be the simplest of elements, its ex- 
pected metallization under pressure fl]] has proved a rather elusive transition. It is now 
accepted that the behaviour of solid and fluid molecular Hydrogen (-H2) may be very 
different at high pressures. Despite considerable experimental efforts with static, room 
temperature compression of solid H2 in diamond anvils beyond pressures of 2 Mbar, 
there is still no compelling evidence for a metallic state Q. The situation is somewhat 
more favourable for fluid H2, since shock compression to 1.4Mbar, and a temperature 
of about 3000K, led to measurements of metallic resistivities ||]. However, theoretical 
interpretation is hampered by the absence of a clear-cut scenario; in particular it is not 
clear whether molecular dissociation precedes ionization or conversely |J. The pres- 
ence of several species, H2, , H, H + and electrons at a "plasma phase transition" || 
complicates a theoretical analysis considerably; to gain a clearer picture of pressure- 
induced ionization, it may be instructive to consider a model system, which would not 
involve molecular dissociation. 

The model system considered in this letter is spin-polarized Hydrogen. If all elec- 
tron spins are assumed to be polarized by a strong external magnetic field B, only 
triplet pair states 3 S can be formed, preventing the binding into H2 molecules. The 
low pressure phase will be made up of H atoms and the only possible scenario upon 
compression will be the ionization of atoms to form an electron - proton plasma, which 
is expected to be crystalline at low temperature and fluid at higher temperatures. A 
rough estimate of the magnetic field needed to spin - polarize the electrons is obtained 
by equating the magnetic coupling energy —fj,sB (where is the magnetic moment 
of an electron) to the difference between the triplet and singlet H-H potential energy 
functions, calculated at the equilibrium distance of the H2 molecule ||; this leads to 
B ~ lO^Tesla. This value exceeds the highest magnetic fields achievable in a laboratory 
by 3 orders of magnitude, but is well within the range of astrophysical situations. The 
present calculation neglects possible orbital effects due to a strong applied magnetic 
field; such effects are expected to be small and their inclusion would lead to a much 
more involved and less transparent calculation. We prefer to think of our model sys- 
tem as a plasma which has been prepared in a spin-polarized state, and is assumed to 
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remain such even when B is switched off. The subsequent calculation will be restricted 
to fluid Hydrogen. 

The thermodynamic properties of the low pressure atomic H j phase may be easily 
calculated from the known triplet pair potential by standard methods of the theory 
of classical fluids. We calculated the atom-atom pair distribution function g(r) from 
the HNC integral equation || , and deduced from it the equation of state via the virial 
and compressibility routes. The resulting excess free energies per atom are plotted 
in Fig. 4 as a function of the usual density parameter r s = a/ao, along the isotherm 
T = W*K; here ao is the Bohr radius, and a = [3/(47r n)] 1 / 3 , where n is the number 
of H atoms per unit volume. There is a thermodynamic inconsistency, typical of HNC 
theory, but the small difference between the "virial" and "compressibility" free energies 
will have no influence on our conclusions. To allow for a meaningful comparison with 
the free energy calculated for the high pressure plasma phase, the free energies shown 
in Fig. 4 contain an electron binding energy contribution of —0.5 a.u. . It is implicitly 
assumed that this binding energy, valid for isolated atoms (i.e. in the limit r s — > oo) 
does not change upon compression up to r s = 2.5, due to overlap and distortion of the 
individual electron Is orbitals. 

A statistical description of the high pressure phase is more challenging. The key 
parameter characterizing the spin-polarized electron component is its Fermi energy 
€f = 2.923/rj? a.u. ;the corresponding Fermi temperature Tp ~ 9.210 5 /r^ K. Along 
the isotherm T = 10 4 -fT considered in the present calculations, the electrons may be 
considered to be completely degenerate (i.e. in their ground state) up to r s m 3. The 
degeneracy temperature of the protons is 2000 times smaller, so that for T = 10 4 -KT, 
the latter may be considered as being essentially classical, down to r s ~ 0.5. The 
proton component is characterized by the Coulomb coupling constant T = e 2 / (aksT) = 
31.56/Vs along the above isotherm, showing that classical Coulomb correlations are 
expected to be strong over the density range 1 < r s < 3 considered in this paper. Note 
that while T decreases as r s increases, the corresponding electron Coulomb coupling 
constant 7 = e 2 /(aep) = 0.342 r s increases. 

In the ultra- high density regime r s < 1, the electron kinetic energy dominates, 
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and the proton and electron components decouple in first approximation ("two-fluid" 
model) ; the weak proton-electron coupling may be treated by linear response theory || , 
suitably adapted to the spin-polarized case. Within linear response, the free energy per 
atom (ion-electron pair) splits into three terms: the ground-state energy of the uniform, 
spin-polarized electron gas ("jellium"), e e , the free energy of protons in a uniform 
neutralizing background (the so-called "one-component plasma" or OCP), focp, and 
the first order correction due to linear screening of the Coulomb interactions by the 
electron gas, Af: 

f = ^f^- = e e + focp + ^f (1) 
where e e (r s ) is taken to be the sum of kinetic (1.754/r s 2 ), exchange (— 0.5772/r s ) 



and correlation [1C] contributions; focp is given by an accurate fit to Monte Carlo 
simulations of the OCP |9|, pi]] ; Af follows from first order thermodynamic perturbation 
theory @: 



A f = ^TZW / S OC p(k)w(k)dk (2) 



1 

2(2^r) 

where Socp(k) is the static structure factor of the OCP (which plays the role of ref- 
erence system). According to linear response theory, w(k) is the difference between 
screened and bare ion-ion pair potentials: 

1 



w(k) 



e(k) 
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(3) 



A- 2 

where e{k) is the dielectric function of the electron gas which we calculated within the 
RPA from the Lindhard susceptibility of a gas of spin-polarized, non-interacting elec- 
trons, supplemented by a local exchange and correlation correction [jH]] . All necessary 
ingredients for the calculation of Af may be found in Q, and the resulting free energy 
curve is shown in Fig. 4. Although linear response cannot, a priori, be expected to be 
quantitatively accurate for r s > 1, it provides a rough estimate of the plasma to atomic 
phase transition, from the intersection of the free energy curves, which is seen to occur 
at r s ~ 1.9. The corresponding transition pressure would be 2.3Mbar. 

However, as r s increases, the ion-electron coupling becomes stronger, and the non- 
linear response of the electron component to the "external" potential field provided 
by the protons, is expected to lower the free energy of the plasma phase. To explore 
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the non-linear regime we have adapted the HNC-DFT formulation of our earlier work 
on (unpolarized) metallic H to the spin-polarized case. Within this formulation, 
proton-proton and proton-electron correlations are treated at the HNC level, which is 
expected to be a good approximation for the long-range Coulomb interactions, while 
the energy of the inhomogeneous electron gas follows from the density functional (E = 
Ne e ): 

E {p{v)\ = E K [p(r)] + E H [p{v)} + E x [p{v)\ + E c [p{v)} (4) 

where p(r) denotes the local electron density, and Ek, Eh, Ex and Eq are the kinetic, 
Hartree, exchange and correlation contributions. For Ek we adopted the Thomas-Fermi 
approximation, corrected by a square gradient term: 

E K W)\ =C K J [p(v)]^dv + ±J l -^fdr (5) 

where Ck = 3(67r 2 ) 2 / 3 /10 a.u., while the choice of 1/9 < A < 1 will be specified below. 
The mean field Hartree term is of the usual form: 

* HWr) ,-!/*/*/M^p (6) 

where Ap(r) = p(r) — n, while: 

E x [p{v)] = C x J[p{v)]^dv (7) 

with Cx = — 3(6/7r) 1//3 /4 a.u. . The correlation contribution Ec\p{v)\ (within the 
LDA) can be found in JlC| ]. This functional yields an explicit form for the electron- 
electron direct correlation function C22(r) (henceforth the indices 1 and 2 will refer to 
protons and electrons respectively) [^2|. The remaining direct and total correlation 
functions cn(r), Ci2(r), hu(r) and h±2{r) are calculated by a numerical resolution of 
the HNC closure equations and the quantum version of the Ornstein-Zernike (OZ) 



relations [ 13 1 , which form a closed set of coupled non-linear integral equations for the 
four functions. 

Solutions were obtained by a standard iterative procedure along the isotherm T = 
W^K and for density parameters in the range 0.5 < r s < 2.5, corresponding to more 
than one- hundred- fold compression of the lowest density state (r s = 2.5), which would 
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correspond to 0.17gr/cm 3 . The temperature is roughly equal to that expected inside 
Saturn, and comparable to temperatures reached in shock compression experiments on 
the NOVA laser facility at Livermore ||. The iterative solutions were first obtained at 
the highest densities (r s = 0.5 and 1), where linear response theory provides reasonably 
accurate intial input. The prefactor A in the square gradient correction to the electron 
kinetic energy functional (0) was adjusted to provide the best match between the HNC- 
DFT result for the local radial density of electrons around a proton, r 2 512 (r) = r 2 [l + 
h\2 an d its linear response prediction, at the highest densities (r s = 0.5 and 1), 
where linear response should be most accurate; g±2 (r) turns out to be rather sensitive 
to A, and the best agreement is achieved for A = 0.18, which is close to the value 1/5 
frequently advocated in electronic structure calculations for atoms [fLOf], 

Results for the local radial density r 2 gvzir) are shown in Fig.l for several values of 
r s . As expected, electrons pile up increasingly at small r as r s increases, and a shoulder 
is seen to develop around r/a ~ 0.4. For comparison, the linear response prediction 
is shown at r s = 1, while at the lowest density (r s = 2.5), an estimate of gviir) in 
the atomic phase is obtained by adding to the electron density in a H atom (namely 
r 2 p(r) = r 2 exp(— 2r)/ir a.u.) the convolution of the latter with the atom-atom pair 
distribution function g(r). HNC-DFT results for the proton-proton pair distribution 
function are shown in Fig. 2 for three densities. As expected, proton-proton correlations 
are seen at first to weaken, as the density decreases, due to enhanced electron screening. 
However at the lowest density (r s = 2.5), weakly damped oscillations build up at long 
distances, which may be indicative of an incipient instability of the proton-electron 
plasma. The atom-atom g(r) at the same density agrees reasonably well with gn(r) 
up to the first peak, but it does not exhibit the long-range correlations in the latter. 
The proton-proton structure factors <Sii(£;) are plotted in Fig. 3 for several values of 
r s . A considerable qualitative change is again seen to occur at the lowest density 
(r s = 2.5), where the main peak is shifted to larger k, while a significant peak builds 
up at k = 0. Such enhanced "small angle scattering" is reminiscent of the behaviour 
observed in simple fluids near a spinodal (subcritical) instability. In fact we were unable 
to obtain convergence of the HNC-DFT integral equations for r s > 2.5, which hints at 
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an instability of the electron-proton plasma at lower densities. This strongly suggests a 
transition to the insulating atomic phase, but the simple density functional used in this 
work cannot properly describe the recombination of protons and electrons into bound 



(atomic) states [14]. 

In order to confirm this scenario, the free energy of the plasma phase should be 
compared to that of the atomic phase. This is easily achieved within the high density 
linear response regime, as shown earlier. However the calculation of the free energy 
in the non- linear regime appropriate for lower densities (r s > 1) is less straightfor- 
ward . In fact the present HNC-DFT formulation provides only one direct link with 
thermodynamics, namely via the compressibility relation ||]: 

lim SWA:) = lim SWA;) = nksT y_ (8) 

fc-»0 fc-»0 

where x T denotes the isothermal compressibility of the plasma. From the calculated 
values of x T , the free energy of the plasma follows by thermodynamic integration, 
starting from a reference state (e.g. r s = 1) for which the linear response estimate is 
expected to be accurate. The resulting "compressibility" free energy curve is plotted 
in Fig. 4. Somewhat unexpectedly it lies above the linear response prediction. An al- 
ternative route to the free energy is via the virial relation for the pressure; only an 
approximate virial expression is known within the present HNC-DFT formulation , 
and the resulting "virial" free energy curve is also shown in Fig. 4. It falls well below the 
"compressibility" free energy curve, thus illustrating the well known thermodynamic 
inconsistency of the HNC closure for Coulombic fluids Any reasonable extrapola- 
tion of the two free energy curves would miss the low density limit —0.5 a.u. by as much 
as twenty per cent. We suggest instead an estimate of the free energy of the plasma 
phase by taking the average of the "compressibility" and "virial" values, despite the 
lack of foundamental a priori justification for doing this. A short extrapolation of the 
resulting curve is likely to intersect or smoothly join on to the free energy of the atomic 
phase just beyond r s = 2.5. An intersection would correspond to a first-order phase 
transition, reminiscent of the "plasma phase transition" of Saumon and Chabrier ||. 
However, due to the uncertainty on the thermodynamics of the plasma phase, a con- 
tinuous transition cannot be ruled out. The transition pressure P would be of the 
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order of 0.5Mbar, well below the current experimental and theoretical estimates for 
the transition of fluid molecular Hydrogen to a conducting state @, || . 

In summary the structure and thermodynamic results derived from an HNC-DFT 
theory of the spin-polarized proton-electron plasma strongly suggest that this plasma 
will recombine into an insulating atomic phase at r s ~ 2.5, for a temperature T = 10 4 X . 
We are presently exploring the behaviour of the system at lower temperatures. 
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Figure 1: The radial electron density around a proton, r 2 512 ( r )? vs r / a ) f° r r s = l(full 
curve), 1.5(dots),2(dash-dotted curve) and 2.5(dashed curve). The triangles are the 
linear response results at r s = 1; the circles correspond to r 2 gn(r) in the atomic 
phase. 
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Figure 2: The proton-proton pair distribution function gu(r) vs. r/a, for r s = l(full 
curve), 2 (dash-dotted curve) and 2.5(circles). The dotted curve is the atom-atom func- 
tion g{r) at r s = 2.5. 
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Figure 3: The proton-proton structure factor Sn(k) vs. ka, for r s = l(full curve), 
1.5 (dash-dotted curve), 2(dots) and 2.5(circles). 




Figure 4: Free energy per proton-electron pair (or atom), minus the ideal proton contri- 
bution, versus the density parameter r s . Full curve: "virial" free energy of the plasma 
phase; dotted curve: "compressibility" free energy for the plasma phase; circles: aver- 
age of the latter two estimates; triangles: linear response result; dash-dotted curves: 
"virial" (upper) and "compressibility" (lower) free energy for the atomic phase. 



